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1 •  Introduct  ion 

The  use  of  uniformly  or  non-unif ormly  distributed 
pseudorandom  numbers  in  systems  simulation,  statistical 
sampling  experiments  and  analytical  Monte  Carlo  work  is  by 
now  well  established.  Numerous  algorithms  sxist  for 
producing  such  numbers  from  various  distributions;  for 
summaries  of  common  techniques,  see  Knuth  [5],  Gaver  and 
Thompson  [2]  or  Ahrens  and  Dieter  [1]. 

The  user  of  pseudorandom  numbers  is  usually  not 
concerned  with  the  details  of  the  algorithm  employed  but 
rather  with  the  results;  a  good  algorithm,  then,  is  one 
which  is  fast,  uses  minimum  computer  memory  and  produces 
numbers  with  satisfactory  statistical  properties.  The 
search  for  statistically  competent  algorithms  for 
pseudorandom  numbers  has  resulted  in  the  specification  of 
many  so-called  "exact"  generators,  that  is  those  whose 
deviation  from  the  true  distribution  concerned  is  the  result 
of  computer  rounding  errors  rather  than  any  defect  in  the 
method  itself.  Such  methods  for  nonuniform  random  numbers 
are  often  based  on  the  assumption  that  "good"  uniform 
numbers  are  available  from  an  independent  generator. 

Exact  generators  for  nonuniform  pseudorandom  numbers 
are  often  quite  complex  and  so  assembly-level  coding  is 
often  resorted  to  when  implementing  thera  in  order  to  meet 
the  computer  time  and  memory  constraints  on  a  good 
algorithm.  An  example  is  the  LLRANDOM  package  developed  at 
the  Naval  Postgraduate  School  by  G.P.  Learraonth  and  P. A.M. 
Lewis  and  described  in  [7];  it  produces  pseudorandom  numbers 


from  uniform,  normal  and  exponential  distributions.  This 
report  describes  an  extension  to  the  LLRANDOM  package  for 
Cauchy  and  gamma  distributed  numbers. 

The  Cauchy_  distribution  has  density  function 


(1) 


f(x)  =  1 


_j.___.  i 


T~+~x* 


-  00  <   X   <   °°    , 


and  distribution  function 


F(x)  =  1   +   1  tan-*  x. 
2 

*•  IT 

While  the  shape  of  the  Cauchy  density  resembles  the  normal 
density,  the  tails  are  much  heavier;  in  fact,  Cauchy  random 
variables  have  no  expectation  and  an  infinite  variance.  The 
density  has  mode  at  zero  and  often  in  applications  the 
variates  are  often  shifted  by  a  location  parameter  T  or 
scaled  by  multiplying  by  a  scale  parameter  S.  Because  of 
the  heavy  tails,  Cauchy  variates  might  find  application  as  a 
"pathological"  case  in  a  systems  simulation  study  as  well  as 
in  statistical  sampling  experiments  for  robust  estimation 
techniques.  See  Chapter  16  of  Johnson  and  Kotz  [4]  for 
further  details  on  the  Cauchy  distribution. 


The  gamma  distribution  with  shape  parameter  A  and  scale 
parameter  s  has  the  density  function 

A  A-1  -sx 
(2)     f(x)  =   s  x   e      , 
— FTAT 

where   r  (A)  is  Euler*s  gamma  function 

A-1  -x 


(3) 


r  (A)  = 


dx 


Note   that    r (n)  =  (n-1) !  when  n  is  a  non-negative  integer 
If  the  random  variable  X  has  density  (2)  then 

E[X]  =  A  /  s   , 


V[X]  =  A  /  s* 

When  A  =  1,  X  has  the  exponential  distribution  while  X, 
suitably  scaled,  has  an  asymptotically  normal  distribution 
as  A  ->  °°  . 

We  note  that  if  X  has  a  r  (A,1)  distribution  then  X/s 
has  a  r  (A,s)  distribution,  so  we  may  set  s  =  1  in  (2)  as 
far  as  the  generating  algorithm  is  concerned.  The  output 
from  the  generator  may  then  be  appropriately  scaled. 

Gamma  random  variables  are  used  in  a  wide  variety  of 
applications:  for  analytical  modeling,  in  reliability 
theory  and  for  statistical  testing  (the  chi-sguared   random 

variable   with   n   degrees   of   freedom   has   the    r(n,   1) 

2   2 

distribution).  See  [6]  or  Chapter  17  of  [4]  for  more 
details. 


II-  U§i=  °f  the  Subroutines 

This  extension  to  LLRANDOM  is  composed  of  two 
independent  IBM  System/360  Assembler-coded  subroutines: 
CAUCHY  for  Cauchy-distributed  variates  and  GAMA  for  gamma 
variates.  The  name  GAMA  was  chosen  so.  as  not  to  conflict 
with  the  IBM  mathematical  library  subprogram  SAMMA  which 
computes  the  gamma  function  (3) . 

The  basic  conventions  for  using  GAMA  and  CAUCHY  are  the 
same  as  in  the  LLRANDOM  package:   the  invoking  statements 

CALL  CAUCHY  (  IX,  X,  N  ) 
and  CALL  GAMA  (  A,  IX,  X,  N  ) 

will  result  in  a  vector  X(1),  ...  ,  X (N)  of  Cauchy  or 
r  (A, 1.0)  pseudorandom  variates,  respectively.  The  arg- 
ument IX  is,  in  both  cases,  an  integer  seed  to  be  used  in 
the  multiplicative  congruential  uniform  generator  employed 
by  LLRANDOM.  IX  should  be  initialized  just  once  in  the 
calling  program  to  some  positive  integer  value  and  should 
not  be  altered  thereafter. 

The  subroutine  GAMA  requires  a  source  for  normal  and 
exponential  deviates;  these  are  obtained  directly  from  the 
LLRANDOM  package  and  so  the  statement  "CALL  OVFLOW"  must 
appear  once  in  the  calling  program  to  initialize  LLRANDOM. 
As  mentioned  previously,  the  output  from  GAMA  must  be  scaled 
if  the  scale  parameter  is  other  than  one;  the  following  set 
of  statements  will  thus  be  required  to  generate  a  vector  of 
100  chi-squared  variates  with  seven  degrees  of  freedom: 

DIMENSION  X  (100) 
CALL  OVFLOW 
IX  =  13726 

•  •  • 

CALL  GAMA  (  3.5,  IX,  X,  100  ) 


DC  50  I  =  1,100 
X(I)  =  2.0  *  X(I) 

50  CONTINUE 

END 

Cauchy  variates  are  also  often  modified  by  location  and 
scale  parameters;  since  no  expectations  exist,  however,  we 
cannot  refer  to  these  parameters  in  terms  of  mean  or 
variance.  Subroutine  CAUCHY  is  completely  independent  of 
LLRANDOM  or  any  other  subroutines  so  that  the  "CALL  OVFLOW" 
statement  is  not  necessary  in  this  case.  To  use  CAUCHY  to 
produce  a  single  variate  C  with  location  parameter  T  and 
scale  parameter  S  we  may  use  the  statements 

•  «  • 

IX  =  217663541 

•  •  • 

CALL  CAUCHY  (  IX,  C,  1  ) 
C  =  S  *  C  +  T 

END 

Just  as  in  LLRANDOM,  linkage  overhead  between  the 
calling  program  and  GAMA  or  CAUCHY  will  be  minimized  if  a 
vector  of  several  variates  is  obtained  at  the  same  time 
instead  of  just  a  single  one.  The  gain  in  this  case  can  be 
as  much  as  50  microseconds  per  variate  in  average  generation 
time,  an  improvement  of  up  to  50X.  In  GAMA,  several 
constants  must  be  calculated  for  each  different  value  of  the 
shape  parameter  A;  these  constants  are  saved  between  calls 
so  that  they  need  not  be  recomputed.  It  will  thus  be  more 
efficient  to  get  several  gamma  variates  with  the  same  shape 
parameter  before  changing  the  A  value,  especially  when  A  > 
3.0   when  the  setup  computations  are   extensive   (see   lines 


174-246  of  the  program  listing)  . 

Note  that  the  techniques  used  in  GAMA  and  CAUCHY  make 
use  of  so-called  rejection  methods  so  that  the  number  of 
uniform  (or  exponential  or  normal)  deviates  needed  to 
generate  a  single  output  deviate  is  random.  When  normal  or 
exponential  deviates  are  required  by  GAMA  from  LLRANDOM  a 
vector  of  1 0  deviates  is  called  for;  since  not  all  of  these 
may  be  used  at  the  time  they  are  generated,  the  balance  are 
saved  for  the  next  call  to  GAMA.  Thus,  reinitializing  the 
seed  IX  to  its  original  value  will  not  in  general  result  in 
an  exact  repetition  of  the  generated  gamma  sequence  since 
the  first  few  deviates  will  use  the  old  normal  or 
exponential  deviates  from  the  previous  sequence.  To  achieve 
an  exact  repetition,  the  generator  must  be  forced  to  repeat 
the  initialization  computations  for  the  desired  A  value;  at 
this  time  any  remaining  variates  from  LLRANDOM  are 
discarded.   An  example  of  this  might  be 

DIMENSION  G  (100) 
CALL  OVFLOW 
IX  =  12345 

•  •  • 

CALL  GAMA  (  A,  IX,  G,  100  ) 

•  •  • 

C      REINITIALIZE  GAMMA  SEQUENCE 
CALL  GAMA  (  1.0,  IX,  G,  1  ) 
IX  =  12345 

CALL  GAMA  (  A,  IX,  G,  100  ) 


END 


CAUCHY  requires  552  bytes  and,  as  mentioned  previously, 
is  completely  independent  of  any  other  subprograms.  CAUCHY 
uses   the   LLRANDOM   multiplicative   congruential   uniform 


generator  but  this  is  coded  in  line  when  needed  so  as  to 
preserve  CAUCHY's  independence.  The  average  generation  time 
per  variate  for  subroutine  CAUCHY  on  a  System/360  Model  67 
under  OS/MVT  was  67.5  microseconds  when  variates  were 
generated  in  vectors  of  100.  The  generation  of  variates  one 
at  a  time  increased  the  average  time  to  119.3  microseconds 
per  variate. 

Subroutine  GAMA  itself  uses  only  1988  bytes  of  memory 
but  since  it  calls  on  LLRANDOM  the  total  core  requirement  is 
9342  bytes: 

GAMA  1988  bytes 

LLRANDOM  6189  bytes 

Required  IBM  Functions  1.1.65  bytes 

Total  9342  bytes 

Timing  the  gamma  generator  on  a  System/360  Model  67  was 
carried  out  using  the  TIME  macro;  Table  1  summarizes  the 
observed  times  as  a  function  of  the  shape  parameter,  A. 
Note  that  since  special  methods  are  employed  when  A  is  0.5, 
1.0,  1.5,  2.0  or  3.0,  the  times  in  these  cases  are 
considerably  shorter  than  times  for  nearby  values  of  A. 


Shape  Parameter 

Algorithm 

Vector  of  100 

Single 

A 

Variates 

Variate 

0.  1 

GS 

324.0 

364.0 

0.3 

GS 

367.0 

402.5 

0.5 

GA 

70.4 

207.7 

0.8 

GS 

439.8 

551.2 

0.9 

GS 

459.0 

611.0 

1.0 

GA 

68.7 

158.9 

1.2 

GF 

300.1 

385.0 

1.4 

GF 

306.  1 

441  .0 

1.5 

GA 

141.7 

215.8 

1.8 

GF 

343.6 

390.8 

2.0 

GA 

142.5 

203.6 

2.1 

GF 

396.1 

450.8 

2.5 

GF 

434.7 

468.5 

2.9 

GF 

444.5 

496.6 

3.0 

GA 

206.7 

237.1 

3.1 

GO 

341.5 

435.8 

3.5 

GO 

336.2 

373.4 

4.0 

GO 

332.4 

420.7 

5.0 

GO 

307.7 

363.2 

8.0 

GO 

293.1 

371.3 

10.0 

GO 

289.4 

312.5 

20.0 

GO 

238.2 

321  .6 

50.0 

GO 

197.7 

284.2 

100.0 

GO 

178.4 

220.0 

1000.0 

GO 

166.7 

177.0 

10000.0 

GO 

136.4 

169.8 

100000.0 

GO 

152.5 

235.8 

Table      1.      Average   generation  times    (microseconds)     for   gamma 
variates    using  subroutine   GAMA. 


III.  Description  of  the  Algorithms 

This  section  describes  the  actual  algorithms  used  in 
CAUCHY  and  GAMA.  An  understanding  of  the  algorithms  is  not 
necessary  for  use  of  the  package  but  they  are  set  forth  here 
both  in  the  interest  of  completeness  and  in  an  effort  to 
document  the  programs  more  fully.  A  single  algorithm 
suffices  for  the  Cauchy  generator  while  GAMA  uses  one  of 
four  algorithms,  depending  on  the  value  of  A. 

In  the  descriptions  which  follow,  the  letters  U,  N  and 
E  (with  or  without  affixes)  represent  uniform,  standard 
normal  and  unit  exponential  pseudorandom  deviates, 
respectively.  The  phrase  "Generate  0"  implies  that  0  is  the 
next  seguential  uniform  variate  in  the  linear  congruential 
sequence;  these  variates  are  generated  as  needed  by  using 
the  same  multiplicative  congruential  scheme  as  used  in 
LLRANDOM.  The  phrases  "Generate  N"  or  "Generate  E"  imply 
that  normal  or  exponential  variates  are  to  be  obtained  by 
linking  directly  to  LLRANDOM. 

A.  Cauchy  Generator 

The  Cauchy  generator  is  a  combination  decomposition- 
rejection  method  (see  Knuth  [5]).  The  Cauchy  density  is 
decomposed,  as  in  Figure   1,   into   three   subdensities:    a 

uniform  density  between  0  and  1  (f  )  ,  a  wedge-shaped  density 

1 

(f  )  and  a  long  tailed  density  (f  ) . 
2  3 


The  uniform  density  f  is  sampled  with  probability 
1/ ^  ;  in  this  case  a  uniform (0,1)  variate  is  returned.  The 
density   f   is  dealt  with  by  using  Marsaglia's  almost-linear 


density  algorithm,  just  as  in  Knuth*s  Algorithm  L  [ 5  ].  The 
density  f  is  sampled  with  probability  1/2  -  1/tt  .  The  tail 
density  f  is  sampled  by  a  rejection  method  with  probability 
1/2.  The  majorizing  density  for  f  is  g(x)  =  1  /  x2 ,  which 
is  the  density  of  the  reciprocal  of  a  uniform  (0,1)  variate. 

Algorithm  C   below   uses  the   fact   that  in  the  prime 

modulus   congruential   random   number   generator    used   in 

LLRANDOH   the  low   order   bits  are  uniformly  distributed  so 

that  b   and  b  select  the  proper  sub-distribution  in  Step  1. 
12 

This   will  not  in  general  be  the  case  for  other  congruential 

pseudo-random  number  generators. 


x), 

i 

1  „ 

2 

\g<x) 

TT     * 

""** f(x) 

1 

f 

2 

TT      " 

1^^ 

fi 

I 

1 
I 

f3 

.5 

1*0 

l'.5 

Figure  1.   Decomposition  of  the  Cauchy  Density  Function. 
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klSL2Li%.k®  0.-        Cauchy   variates. 

1.  (Select      subdensity)         Generate   U,    setting    aside   the   two 

low    order   bits   b      and  b    .      If   b      =    1,    go   to   step   6. 

1       2        1 

2.  (Sample   box)      If  U    <    0.6366197724  =  ,2/tt  ,  generate  a 

*  * 

new  variate  U  ,  set  x  =  U   and  go  to  Step  8. 

3.  (Sample  wedge)    Generate  new  variates  U   and  U    .      If  0 

1       2        1 

>  U  ,  exchange  U   and  U  .   Set  x  =  U  . 
2  12  1 

4.  (Easy   rejection)     If  U   <  0.8284271247  =  2  /I  -    2,    go 

2 

to    Step    8. 


5.  (Hard   rejection)         If    U      -    U      <    1         x*     (    2/2   -    2) ,    go   to 

2     1    T-~x* 

Step  8,    otherwise  go  back  to  Step  3. 

6.  (Sample  tail)    Set  x  =  1  /  U. 

*        * 

7.  (Tail  rejection)    Generate  a  new  variate  U  .   If   U    < 


x£     go   to  Step  8,  otherwise  generate  a  new  U  and  go 
T  +  x"2 

back  to  Step  6. 

8.   (Random  sign)    If  b  =1  set  x  =  -  x.   Deliver  x  as  the 

2 

generated  deviate. 


It  should  be  noted  that  there  are  several  other  methods 
for  generating  Cauchy  variates:  the  ratio  of  independent 
standard  normal  deviates  has  the  Cauchy  distribution,  as 
does  the  quantity 

X  =  tan  [  *   (u  -  1)  ]r 

7 

where  U  is  unif orm  (0, 1) .   These  methods  are   both   substant- 
ially slower  than  algorithm  C,    but  another  new  method  has  an 


11 


average  time  comparable  to  Algorithm  C  and  is  much  easier  to 
program.  This  second  method  requires  an  average  of  2.55 
uniform  random  variates  per  Cauchy  variate  (as  compared  with 
2.47  for  algorithm  C)  and  it  needs  about  69  microseconds  per 
variate  on  the  System/360  Model  67.  It  is  possible, 
however,  that  Algorithm  CR  will  be  better  than  algorithm  C 
in  some  other  implementation. 

The  method  is  essentially  the  technique  devised  by  von 
Neumann  to  generate  a  random  variate  sin  U,  where  U  is 
uniform  between  0  and  2tt  .  Such  variates  are  used  in  the 
polar  method  for  generating  normal  random  variables  [8].  It 
does  not  seem  to  have  been  recognized  that  the  method  also 
generates  tan  U,  which  is  the  required  Cauchy  variate. 

Algorithm  CR.   Cauchy  variates,  ratio  method. 

1.  (Get  uniforms)    Generate  U   and  U  .  Set  Y  =  2  U      1 

12       11 

and  Y   =  2  U   -  1. 

2      2 

2.  (Rejection  test)    If  Y*  +  Y2  >  1  go  back  to  Step  1. 

12 

3.  (Take  ratio)         Deliver   x   =   Y      /  Y    . 

12 


B.   Gamma  Generator  GS:  A  <  1.0 

—————   — —  —  —  — — 7r  —  — '   — — —   —  —   — — — 

This  method  is  due  to  Ahrens  and  is  set  forth  in  [ 1 ]. 
It  is  applicable  only  to  values  of  A  less  than  one  and  is 
markedly  superior  in  execution  time  to  the  method  of  Johnk 
[3],  which  is  the  usual  technique  for  generating  variates  of 
this  type. 

The  method  is  a  rejection  method  employing  two 
different  tests,  one  of  which  is  chosen  at  random  for  any 
given   variate;    the   power   transform  of   a   uniform  (0,1) 


12 


VA 

vanate,   U    ,   is   tested  in  the  region  0  <  x  <  1,  while  a 

suitable  exponential,  E,  is  tested  when  x  >  1.  The 
advantage  of  this  method  lies  in  the  limited  use  of  the 
library  subprograms  for  the  exponential  and  logarithm; 
average  times  range  from  300  to  400  microseconds  as  compared 
with  600  to  800  for  Johnk's  method.  Further  discussion  and 
proofs  may  be  found  in  [1]. 


Algorithm  GS.    Gamma  variates,  A  <  1.0. 

1.  (Select   rejection  test)    Generate  U  and  generate  E  and 

set  P  =  e  +  A  U.   (Note  that  "e"   is   the   base   of   the 

e  "" 

natural   logarithms.)    If  P  <  1  go  to  Step  2,  otherwise 

go  to  Step  3. 

VA 

2.  (Small   x   test)      Set  x  =  P    .   If  x  <  E,  deliver  x, 

otherwise  go  back  to  Step  1. 

3.  (Large   x   test)    Set  x  =  -  In  [  1  {  e  +  A  -  P  }  }.   If 

T.  e 

(1  -  A)  In  x  <  E,  deliver  x,  otherwise  go  back   to   Step 
1. 

C.   Gajnma  Generator  GF:  1 . 0  <  A  <  3.0 

A  thus-far  unpublished  method  devised  by  Professor  G.S. 
Fishman  of  North  Carolina  University  was  communicated  to  the 
authors  in  private  correspondence.  It  is  valid  for  any  A  > 
1.0  but  its  efficiency  in  terms  of  average  time  goes  down  as 

v/A"   so   it   is  applied  in  GAMA  only  in  the  range  where  it  is 

superior  to  the  Dieter-Ahrens   method   GO   described   below. 
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The   method   is  a   rejection   method  based  on  the  following 
theorem. 

Theorem   Let  U  be  a  uniform  (0,1)  random  variable  and  let  E 
be  an  exponential  random  variable  with  mean  A.   Let 

A-1   -x(1-1/A)  -  (A-1) 
g(x)  =  c  x  3     e 

If  g (E)  >  U,  then  E  has  conditionally  the  gamma  distribution 
with  shape  parameter  A,  i.e. 

A- 1  -x 
f  (  x  |  0  <  g  (E)  )  =  x   e 

E*  —  ?  in- 

Proof: 


Unconditionally,      E      has      density      h  (x) 


1      e 

A* 


x/A 


Therefore, 


(4)  fE<x|U<g(E)    )    =        hlxJ.^P£lU<a(EHE=xL       . 

Now   since  U    is  uniformly   distributed, 

Pr{U<g  (E)  |E=x}    =   g(x) 

as   long   as  0  <  g  (x)  <  1;  that  this  is  true  for  every  x  >  0 
may  be  readily  verified  by  elementary  calculus.   Therefore, 


(5) 


Pr{U<g(E)}    =    E[    Pr{U<g(E)|E}     ] 
-  g  (x)     h  (x)     dx 

Jo 

A-1    -A 
=      r  (A)     e         A 

=    C(A) 


Thus,    in    view   of    (4)  , 
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f_(x|U<g(E)  )  =  hix^  q[x] 


E 


A-1  -x 
=  x   e 


The  efficiency  of  the  generator  is  governed  by  the 
probability  that  a  given  variate  will  pass  the  rejection 
test,  U  <  g (E) ;  from  (5)  it  will  be  seen  that  this 
probability   is   just   C(A)  .    When   A  is  large  we  have  from 

■ 

Stirling's  approximation  that  C(A)  =   /"Zjr  ,      so   that   the 

>/  A*  e7 

method   becomes  more  inefficient  with  increasing  A,  as  noted 

above. 

A  slight  modification  to  the  method  suggested  by  the 
theorem  improves  the  efficiency  slightly  and  we  obtain 

Ai92£i:ill§  £F.   Gamma  variates,  1.0  <  A  <  3.0. 

1.  (Generate  exponentials)      Generate   two    independent 

exponential  variates,  E   and  E  . 

1  2 

2.  (Rejection  test)    If  E   <  (A-1)  (E   -  In  E   -   1)   then 

2  11 

go  back  to  Step  1. 

3.  (Acceptance)    Deliver  x  =  A  E  . 

D.   Gamma  Generator  GOz    A  >  jkO 

This  method  was  originally  developed  by  Dieter  and 
Ahrens  and  is  fully  described  in  [  1 ]  together  with  several 
other  gamma  generation  techniques.   Algorithm   GO   does  not 
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suffer  the  usual  drawback  of  growing  less  efficient  in 
generation  time  with  increasing  A;  in  fact,  the  method  is 
more  efficient  for  larger  A  values. 

The  basi-c  idea  here  is  to  take  advantage  of  the 
asymptotic  normality  of  the  gamma  distribution  by  doing  most 
of  the  sampling  from  a  normal  distribution;  the  right  hand 
tail  is  sampled,  when  necessary,  using  a  rejection  method 
with  the  exponential  distribution.  The  method  can  be 
applied  to  values  of  A  greater  than  2.533,  but  it  is  not  as 
efficient  as  Fishman's  technigue  for  A  <  3.0. 

As  mentioned  previously,  this  algorithm  requires  the 
computation  of  several  constants  which  depend  only  on  A  and 
which  may  be  saved  between  calls;  these  calculations  are 
described  in  step  0  of  the  specification  below.  Further 
discussion,  illustrations  and  proofs  are  given  in  [ 1  ];  the 
version  of  GO  here  differs  in  a  few  minor  details  from  the 
original  Dieter  and  Ahrens  technique. 

Algorithm  GO.   Gamma  variates,  a  >  3.0. 

0.  (Calculate  constants)    Compute: 

m  =  A  -  1 ; 

s2  =   /"8I  +  A;        s  =  /s*; 

d  =  /Ss?;  b  =  d  ♦  m; 

w  =  s2  /  m  -  1;        v  =  2s2  /  (m  /I   )  ; 

c  =  b  ♦  In  s  d  -  2m  -  3.7203285. 

1.  (Select   normal/exponential)  Generate  U.  If   U  < 
0.0095722652  go  to  Step  7. 

2.  (Normal  sampling)    Generate  N  and  set  x  =  sN  +  m. 

3.  (Check  trial  value)  If  x  <  0  or  x  >  b  go  back  to  Step  2, 
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otherwise  generate   a  new  variate  U  and  set  S  =  N2  /  2. 
If  N  >  0  go  to  Step  5. 

4.  (Left-hand  rejection)    If  0  <  1  ♦  S  (vN  -  w)  go  to  Step 
9,  otherwise  go  to  Step  6. 

5.  (Right-hand  rejection)    If  U  <  1  -  wS  go  to  Step  9. 

6.  (Final  normal  rejection)    If  In  U  <  m  In  x  +  m  -  x  +   S 

m 

go  to  Step  9;  otherwise  go  back  to  step  1. 

7.  (Exponential)   Generate  E   and  E   and  set  x  =  b(1  +  E  /d)  . 

1       2  1 

8.  (Exponential   rejection)     If  m  (x  -  In  x  )  +  c  >  E   go 

B      m  2 

back  to  Step  1 . 

9.  (End)    Deliver  x  as  the  gamma  variate. 


E.   Ad  Hoc  Gamma  Generators 

This   set  of  algorithms  is  based  on  the  well-known  fact 

that  the   sura   of   independent   gamma   variates   with   shape 

parameters   A    and   A    and   equal  scale  parameters  has  the 
1         2 

gamma  distribution  with  shape  parameter   A  +   A    and   scale 

1     2 

parameter   equal   to   that   of   the   summands.    We  may  thus 

generate  a  gamma  variate  with  integer  shape  parameter  K  by 
taking  the  sum  of  K  independent  exponentials.  This  will  be 
more  efficient  than  the  previously  discussed  methods 
(Algorithms  GF  and  GO)  for  moderate  values  of  K;  for  the 
System/360  we  take  K  <  3  to  apply  this  ad  hoc  technique. 

An  obvious  extension  to  this  method  is  to  allow  for 
half-integral  values  of  A  by  making  use  of  the  fact  that  the 
square  of  a  standard  normal  random  variable  has  the 
chi-squared  distribution  with  one  degree  of  freedom,  i.e. 
N2/2  has  the  gamma  distribution  with  unit  scale  parameter 
and  A  =  0.5.   We  use  this  extension  for  A  =  0.5  or  1.5. 
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The  resulting  algorithm  is  then 

Alg,orithin  GA.  Gamma  variates,  integral  or  half-integral 
shape  parameter  A. 

1.  (Find  K)  Set  K  =  [A]f  where  [A]  denotes  the  integral 
part  of  A.  Set  X  =  0.  If  A  -  K  =  0.5  set  L  =  1;  if  A  - 
K  =  0.0  set  L  =  0;  otherwise  Stop.  (If  the  algorithm 
stops,  an  incorrect  A  value  has  been  used.) 

2.  (Generate  exponentials)     If  K   =   0   go  to   Step   3, 

otherwise  generate   K   exponentials  E  ,  ...  ,  E   and  set 

1         K 

X  *  E   +  .  . .  ♦  E  . 
1  K 

3.  (Generate  normal)  If  L  =  0  go  to  Step  4  otherwise 
generate  N  and  set  X  =  X  +  N2/2. 

4.  (Deliver  X)    X  is  the  desired  variate. 
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IV.   Summary  and  Comments 


This  work  provides  a  convenient  and  useful  extension  to 
the  LLRANDOM  package,  especially  for  users  interested  in 
statistical  and  reliability  theory  applications  of  digital 
simulation.  The  combination  of  the  most  efficient  known 
gamma  generation  technigues  with  the  new  Cauchy  method  gives 
exceptionally  good  time  characteristics  at  some  cost  in 
computer  memory  utilization. 

The  work  may  be  extended  at  once  to  the   generation   of 

several   other   types  of  random  variables.   For  example,  the 

beta  distribution  with  parameters  A  and  B  may  be  sampled   by 

talcing   gamma  variates   X    and   X    with   respective  shape 

1         2 

parameters  A  and  B  and  delivering 

•   Z  =  X   /  (  X   +  X   ) 
1       1     2 

as  a  beta  variate.   In  this  case   considerable   overhead   in 

GAMA   can   result  from  shifting  the  shape  parameter  back  and 

forth  between  A  and  B;  for  this  reason  obtaining  vectors   of 

gamma  variates  X   and  X   is  recommended,  as  in  the  following 

12 

example: 

DIMENSION  X1(50),  X2  (50)  ,  Z(50) 

CALL  GAMA  (  A,  IX,  X1,  50  ) 
CALL  GAMA  (  B,  IX,  X2,  50  ) 
DO  405  I  =  1,50 

Z(I)  =  X1  (I)  /  (  X1(I)  +  X2(I)  ) 
405  CONTINUE 


END 
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The  t-Distribution  may  be  sampled  as  the  ratio  of  a 
standard  normal  and  an  independent  chi-squared  random 
variate,  while  the  F~ Distribution  may  be  obtained  by  taking 
the  ratio  of  two  independent  chi-squared  variates  divided  by 
their  respective  degrees  of  freedom.  (See  pages  4  and  5  for 
an  example  of  the  generation  of  chi-squared  variates.) 
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